188_CRISPR off-target分析¶

# 安装Cas-OFFinder
# 下载：http://www.rgenome.net/cas-offinder/

# 准备输入文件
cat << 'EOF' > input.txt
/path/to/hg38.fa
NNNNNNNNNNNNNNNNNNNNNGG
ACGTACGTACGTACGTACGTNNN 5
EOF
# 第1行：参考基因组路径
# 第2行：PAM模式（NGG for SpCas9）
# 第3行：gRNA序列 + 最大错配数

# 运行Cas-OFFinder
cas-offinder input.txt C output.txt
# C = CPU模式, G = GPU模式

# 输出格式
# Bulge_type  crRNA  DNA  Chromosome  Position  Direction  Mismatches  Bulge_Size
head output.txt

# 高级用法：允许DNA/RNA bulge
cat << 'EOF' > input_bulge.txt
/path/to/hg38.fa
NNNNNNNNNNNNNNNNNNNNNGG
ACGTACGTACGTACGTACGTNNN 4 1
EOF
# 最后两个数字：最大错配数 最大bulge大小

# CRISPOR命令行版
# 安装
pip install crispor

# 运行（也可使用Web版：http://crispor.tefor.net/）
crispor.py \
    hg38 \
    ACGTACGTACGTACGTACGT \
    NGG \
    --output crispor_output.tsv

# 输出包含：
# - MIT特异性分数（0-100，越高越特异）
# - CFD分数（每个脱靶位点的切割频率预测）
# - 脱靶位点列表及详细信息
# - gRNA效率评分

# Python中使用CRISPOR API
import requests

def query_crispor(sequence, genome='hg38', pam='NGG'):
    """查询CRISPOR网站获取gRNA评估"""
    url = "http://crispor.tefor.net/crispor.py"
    params = {
        'seq': sequence,
        'org': genome,
        'pam': pam
    }
    # 注意：实际使用应遵循CRISPOR的API使用条款
    response = requests.get(url, params=params)
    return response.text

# 解析Cas-OFFinder结果
import pandas as pd

offtargets = pd.read_csv("output.txt", sep="\t", header=None,
                          names=["BulgeType", "crRNA", "DNA", "Chrom",
                                 "Position", "Strand", "Mismatches", "BulgeSize"])

# 统计脱靶分布
print("按错配数分布:")
print(offtargets['Mismatches'].value_counts().sort_index())

# 过滤高风险脱靶（≤3错配）
high_risk = offtargets[offtargets['Mismatches'] <= 3]
print(f"\n高风险脱靶位点数（≤3错配）: {len(high_risk)}")
print(high_risk[['Chrom', 'Position', 'Strand', 'Mismatches', 'DNA']])

# CFD分数计算
import numpy as np

# Doench 2016论文中的CFD权重矩阵（简化版）
# 每个位置每种错配的权重
def calc_cfd_score(guide, offtarget):
    """计算CFD分数"""
    assert len(guide) == len(offtarget) == 20  # 不含PAM

    # 位置权重（近PAM端权重更高）
    # 实际CFD矩阵更复杂，这里简化
    position_weights = {
        1: 0.0, 2: 0.0, 3: 0.014, 4: 0.0, 5: 0.0,
        6: 0.395, 7: 0.317, 8: 0.0, 9: 0.389, 10: 0.079,
        11: 0.445, 12: 0.508, 13: 0.613, 14: 0.851, 15: 0.732,
        16: 0.828, 17: 0.615, 18: 0.804, 19: 0.685, 20: 0.583
    }

    score = 1.0
    for i in range(20):
        if guide[i] != offtarget[i]:
            pos = i + 1
            # 错配降低分数
            score *= (1 - position_weights.get(pos, 0.5))

    return score

# 示例
guide = "ACGTACGTACGTACGTACGT"
ot1 = "ACGTACGTACGTACGTACGA"  # 1个错配在位置20
ot2 = "ACGTACGTACGTACGTAAGT"  # 1个错配在位置17

print(f"OT1 CFD score: {calc_cfd_score(guide, ot1):.4f}")
print(f"OT2 CFD score: {calc_cfd_score(guide, ot2):.4f}")

# 主要实验检测方法对比

## GUIDE-seq（Genome-wide Unbiased Identification of DSBs Evaluated by Sequencing）
- 原理：将短双链寡核苷酸（dsODN）整合到Cas9切割位点
- 优势：全基因组无偏检测，灵敏度高
- 限制：需要在活细胞中进行

## CIRCLE-seq（Circularization for In vitro Reporting of Cleavage Effects by Sequencing）
- 原理：体外环化基因组DNA后用Cas9切割
- 优势：无需细胞转染，灵敏度极高
- 限制：体外条件可能高估脱靶

## Digenome-seq（Digitized genome sequencing）
- 原理：体外用Cas9消化基因组DNA后测序
- 优势：简单直接
- 限制：需要高深度测序

## DISCOVER-seq
- 原理：检测MRE11蛋白在切割位点的结合
- 优势：检测细胞内真实切割
- 限制：需要特殊抗体和ChIP-seq

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np

# 可视化gRNA与脱靶位点的比对
def visualize_alignment(guide, offtargets_df, output_file="offtarget_vis.png"):
    fig, ax = plt.subplots(figsize=(14, len(offtargets_df)*0.5 + 2))

    guide_display = list(guide)
    y_pos = len(offtargets_df)

    # 画gRNA
    for i, nt in enumerate(guide_display):
        ax.text(i, y_pos + 1, nt, ha='center', va='center', fontsize=10,
                fontweight='bold', color='black',
                bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
    ax.text(-2, y_pos + 1, 'gRNA', ha='right', va='center', fontweight='bold')

    # 画每个脱靶位点
    for idx, (_, row) in enumerate(offtargets_df.iterrows()):
        ot_seq = list(row['DNA'][:20])  # 不含PAM
        y = y_pos - idx

        for i, (g, o) in enumerate(zip(guide_display, ot_seq)):
            color = 'lightgreen' if g == o else 'salmon'
            ax.text(i, y, o, ha='center', va='center', fontsize=9,
                    bbox=dict(boxstyle='round', facecolor=color, alpha=0.7))

        label = f"OT{idx+1} ({row['Mismatches']}mm) {row['Chrom']}:{row['Position']}"
        ax.text(-2, y, label, ha='right', va='center', fontsize=8)

    ax.set_xlim(-8, 22)
    ax.set_ylim(-0.5, y_pos + 2.5)
    ax.axis('off')
    ax.set_title('Off-target Alignment Visualization', fontsize=14)

    plt.tight_layout()
    plt.savefig(output_file, dpi=300, bbox_inches='tight')

# 使用示例
# visualize_alignment("ACGTACGTACGTACGTACGT", high_risk)

# Cas-OFFinder搜索脱靶（GPU加速）
cas-offinder input.txt G output.txt

# CRISPOR Web API（推荐Web版）
# http://crispor.tefor.net/

# FlashFry（快速gRNA设计+脱靶预测）
java -jar FlashFry.jar index --tmpLocation tmp/ --database hg38.fa --output hg38_db
java -jar FlashFry.jar discover --database hg38_db --fasta target.fa --output guides.tsv
java -jar FlashFry.jar score --input guides.tsv --output scored.tsv --scoringMetrics doench2014,doench2016,dangerous

# CRISPRitz（允许bulge的脱靶搜索）
crispritz.py search hg38/ guide.txt PAM.txt output/ --mm 6 --bDNA 2 --bRNA 2 -t